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@ Method for determining petrophysical properties of a subterranean layer. 

@ A method is provided for determining petrophysical properties associated with various lateral locations of a 
subterranean layer. The method employs velocity and density log data, corresponding to a preselected 
reference lateral location, in combination with seismic data to determine a range of values of at least one 
petrophysical property associated with a desired lateral location ("nonreference") offset from the reference 
lateral location. 
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Background of the Invention 

This invention relates to a nnethod for deternnining petrophyslcal properties of a subten'anean layer 
which employs both seismic data and log data. 

5 It is well known to employ logs, such as wireline well logs, to determine valuable petrophyslcal 
properties associated with a subterranean layer of interest. Petrophyslcal properties, such as porosity, water 
or hydrocarbon saturation and llthological composition (i.e. shale or sand), provide valuable information in 
determining the presence and extent of hydrocarbons in the layer of interest. However, such logs are very 
limited in areal extent to only about 6-12 inches around a borehole in which measurements are tal<en. 

10 Moreover, obtaining logs such as wireline well logs can be time consuming and expensive In requiring 
drilling of a borehole for each such well log. 

Petrophyslcal properties of a layer of interest can vary widely at different locations. Therefore, accurate 
determination of variations in properties over a large area of a layer are not practical by use of wireline well 
logs, since such determination would require many, possibly hundreds or thousands, of such well logs. 

75 Seismic prospecting is effective in estimating depths to subterranean layers, and is cost effective in 
surveying a large area, but the resulting seismic data provides insufficient information to make accurate 
detenminations of. for example, the extent and amount of hydrocarbons in a hydrocarbon-containing 
subterranean layer (hydrocarbon reservoir). 

20 Summary of the invention 

It is, therefore, an object of the invention to provide a cost effective method capable of determining 
petrophyslcal properties associated with a subterranean layer of interest at any lateral location thereof. 

The above object is realized by a method described herein as a series of ten steps employing both 
25 seismic and log data. The method permits accurate determination of petrophyslcal properties of a layer of 
interest at any desired lateral location thereof. The method Is cost effective insofar as it requires minimal log 
data (as little as a single set of data corresponding to a single lateral location) to be employed with the 
seismic data. 

30 Detailed Description of the Drawings 

FIG. 1 is a schematic illustration of a cross-section of the earth which shows a subterranean layer of 
interest having a hydrocarbon reservoir extending therethrough. This FIGURE also schematically Illustrates 
seismic and well log equipment for collecting data employed in the method of the invention. 
35 FIG. 2 is a depiction of a portion of a seismic section, where such seismic section is composed of a 
plurality of seismic traces respectively corresponding to various lateral locations of the layer of interest in 
FIG. 1. 

FIG. 3 schematically illustrates the manner in which a plurality of logs, corresponding to a particular 
lateral location of the layer of interest, are employed in production of a synthetic seismogram. 
40 FIG. 4A shows logs which correspond to a reference lateral location in FIG. 1, and FIGS. 4B and 4C 
show modified versions of the logs in FIG. 4A which are employed in the method of the invention. 

D tailed Description of the invention 

45 The method of the invention will now be described in detail in terms of a simple embodiment employing 
a single reference lateral location and a single reference log-par. It should be understood, however, that the 
method could employ multiple reference lateral locations and associated reference log-pairs, as is 
employed in a subsequent example described herein. The number of reference lateral locations and 
reference log-pairs employed depends upon the size of the area for which petrophyslcal properties are 
50 desired, and the extent to which the geology of such area varies. 

Various terms as used herein and in the appended claims are defined as follows. 
A "lateral location", as such term is used herein, is defined by a vertical line wherein different lateral 
locations are horizontally spaced from one another. 

The term "log" refers to a set of at least one data point or series of data points, expressible in terms of 
55 a curve or function of depth or time, representative of a particular physical parameter associated with a 
subterranean formation and obtained by any means, unless a particular means is specified. 

The term "petrophyslcal property" means any property of a subten'anean layer which is related to the 
presence (or lack thereof) and/or amount of hydrocarbons in such layer. Examples of petrophyslcal 
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properties include, but are not limited to, layer thickness, porosity, lithological composition (i.e. volume 
percent of shale or sand), water or hydrocarbon saturation, and any properties derived from one or more of 
the aforementioned properties such as gross pay thickness, net pay thickness, net pore thickness, 
hydrocarbon pore volume, and net porosity. The latter-mentioned five properties are discussed and defined 
5 below with reference to the FIGURES. 

The term "pay" as used in some of the above terms is a shorthand term for hydrocarbon-bearing or 
containing in what is deemed to be sufficient quantities at a given depth to justify commercial production. 

The term "reflection coefficient" Is a parameter which Is calculated from velocity and density as is 
further discussed below, and such term as used by itself refers to a set of at least one value or a set of a 
70 series of values which can be expressed in terms of a curve or function of time or depth. In the latter case 
of a set of a series of values, the term "reflection coefficient series" will sometimes be used. 

The term "crosscorrelation" is well known to those skilled In the art, and refers to a technique of 
measuring the similarity of two waveforms. When normalized, a crosscon^elatlon value of one Indicates a 
perfect match and a value near zero indicates little correlation. 
15 The term "difference mismatch error" as such term is applied to the comparison of a first waveform 
and a second waveform, where each waveform has corresponding data points at predetermined sample 
intervals, is the sum of the absolute values of the differences in amplitudes between corresponding data 
points of the first and second waveforms divided by the number of samples within a designated comparison 
window. 

20 The various steps of the method, which can be performed in a different order if desired, will now be 
described with respect to determination of petrophysical properties of a hydrocarbon reservoir schematically 
illustrated in FIG. 1. There is shown in FIG. 1 a cross-section of the earth at 10 having multiple subterranean 
layers. One of such layers shown at 12 has a primarily hydrocarbon-bearing section 12a, designated as the 
hydrocarbon reservoir, and primarily water wet sections 12b and 12c. The lateral extent of reservoir 12a 

25 ranges from a first lateral location at LI to a fifty-fifth lateral location at L55. Every fifth lateral location is 
shown, with the exception of the twenty-eighth lateral location at L28, which is the reference lateral location 
sometimes denoted hereafter as "reference L28" or simply the "reference lateral location". The other lateral 
locations are nonreference lateral locations sometimes denoted hereafter by the phrase constituted by 
"nonreference" followed by a particular location number or by simply "nonreference lateral locatlon(s)", 

30 . 

1. Provide Reference Seismic Trace and Nonreference Seismic Traces 

The reference seismic trace, designated as T28 in FIG. 2 and sometimes referred to hereafter as 
"reference T28" or simply the "reference seismic trace", corresponds to the reference lateral location, 

35 reference L28. The reference seismic trace, reference T28, can be conventionally obtained by generation of 
at least one seismic pulse at the surface, which travels down to the boundaries of reservoir layer 1 2 so as 
to be reflected by such boundaries and received/detected by one or more receivers. Such a surface 
seismic arrangement is schematically illustrated in FIG. 1. A seismic source is indicated at 14 and a seismic 
receiver at 16. The midpoint between such a source and receiver is at reference L28 such that rays 

40 associated with the seismic pulse, indicated at 18 and 20, are reflected by the upper and lower boundaries 
of layer 12 at reference L28 and accordingly received and detected by receiver 16 to produce reference 
T28. Reference T28 includes the two reflection events corresponding to the respective upper and lower 
boundaries of layer 1 2. 

Of course, in actual practice, a plurality of source-receiver pairs, having a common midpoint at 
45 reference L28, would be employed to obtain a plurality of seismic traces, which would then be corrected for 
normal moveout and stacked to obtain a single composite trace such as reference T28 shown In FIG. 2. 

Similarly, each of the nonreference seismic traces, nonreference T1-T27 and T29-T55, are obtained by 
employing source-receiver pairs not shown. Nonreference T1-T27 and T29-T55 respectively correspond to 
reference L1-L27 and L29-L55, and each Include a pair of reflection events respectively con'esponding to 
50 the upper and lower boundaries of layer 12. 

2. Provide Reference Log-Pair 

A reference log-pair, comprising a velocity log and a density log, is provided for layer 12 at reference 
55 L28. Such reference log-pair should also have associated therewith at least one known petrophysical 
property which is desired to be determined for the various nonreference lateral locations. Velocity can be 
expressed for the velocity log in, for example, feet/second or in the reciprocal form of microseconds^oot 
(called a "sonic log"). 
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The reference log-pair can, according to certain broad aspects of the invention, be obtained by any 
means including, but not limited to, wireline well logs, other types of well logs such as those obtained from 
surface analyses of core samples or cuttings, vertical seismic profiles, and any combinations or derivations 
thereof. It is preferred, however, to obtain the reference log-pair in the manner described below, employing 

5 wireline well logs and derivations of such logs. In FIG. 1 , a well logging tool is schematically indicated at 22 
as being suspended in a borehole 24 by a suitable wireline so as to be positioned at reference L28 between 
the upper and lower boundaries of layer 12. An appropriate type of tool, depending on the parameter being 
measured, is employed to make measurements through layer 12 at various depths at predetermined 
increments. Such increments can range from a few inches to several feet, but are typically about i foot. 

10 According to a preferred embodiment, velocity (represented by "v", i.e. ft./sec.) and density (repre- 
sented by "p", i.e. g/cm^) wireline well logs are taken for layer 12 between the upper and lower boundaries 
of such layer. Such wireline well logs are employed In combination with other appropriate wireline well logs, 
such as gamma, spontaneous potential, and resistivity well logs, to derive a porosity (represented by 
the fraction, i.e. in %, of the total volume of formation material which is pore volume) log in a manner well 

15 known to those skilled in the art. A saturation (i.e. "Sw", the fraction. I.e. in %, of the pore volume of the 
fomrtation material which Is occupied by water) log is similarly derived from resistivity or induction wireline 
well logs, and a lithological composition log (i.e. "Vsh", the fraction, i.e. In %, of the total volume of 
formation material which is shale) is derived from a gamma or spontaneous potential wireline well log. Such 
Vsh, and Sw logs are shown at 26, 28, and 30, respectively, in FIG. 3. Each of Vgh, and S^ are 

20 petrophyslcal properties associated with the reference lateral location, reference L28, from which other very 
useful petrophyslcal properties can be derived as will be explained in another step. 

New velocity (v) and density (p) logs, such as those shown in FIG. 3 at 32 and 34, respectively, are 
derived from the Vsh, ^, and Sw logs in a manner described in detail in a subsequent example, thereby 
providing the desired reference log-pair. 

25 

3. Determine Reference Reflection Coefficient 

From the velocity and density logs 32 and 34 of the reference log-pair, the reference reflection 
coefficient is preferably determined as a series of values, schematically indicated at 36 In FIG. 3. from the 
30 well known formula 



35 

wherein a reflection coefficient value is calculated from such formula for each corresponding data point-pair 
at predetermined sample intervals. A data point-pair from velocity log 32 would comprise V2 and vi 
corresponding to data points separated by the predetemnined sample Interval, and a con^esponding data 
40 point-pair from density log 34 would comprise p2 and pi con'esponding to data points separated by the 
predetemiined sample interval. 

4. Provide Seismic Wavelet 

45 At least one seismic wavelet Is provided which is representative of the seismic pulse(s) at the layer of 
interest, layer 1 2, and which when convolved with the reference reflection coefficient series determined in 
step 3 produces a reference synthetic seismogram which approximates the reference seismic trace, 
reference T28. 

A seismic wavelet can be derived by various techniques. For example, the actual seismic pulse(s) can 
so be measured at the surface and then corrected by appropriate data processing to account for distorting 
effects of the earth between the surface and layer 12. Or, more preferably, the wavelet is extracted from a 
line of seismic traces (I.e. seismic section portion) con'esponding to lateral locations closely adjacent to 
reference 1-28. such as T1-T55. 

The seismic wavelet produced as discussed above is convolved with the reference reflection coefficient 
55 to produce a reference synthetic seismogram. If comparison (over corresponding comparison windows 
including reflection events corresponding to the upper and lower boundaries of layer 12) of the thus 
produced reference synthetic seismogram with the reference seismic ti^ace, T28. by any suitable compari- 
son technique (such as crossconrelation), passes at least one predetermined matching threshold (i.e. a 
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minimum crosscorrelation value of, for example, 0.90 or 0.95), such seismic wavelet can be employed in 
subsequent step (7). Such comparison will necessitate approximate alignment of the reference seismic 
trace with the reference synthetic seismogram, and may also preferably involve shifting of either with 
respect to the other to determine the alignment which gives the best match. If such comparison does not 

5 pass the predetermined matching threshold, the seismic wavelet can be modified in shape by modifying 
wavelet frequency, amplitude, and/or phase until such matching threshold is passed. Alternatively, a 
predetermined number of wavelets can be generated having slightly different shapes, and the wavelet 
selected for use in step (7) which, when convolved with the reference reflection coefficient series, produces 
the reference synthetic seismogram which passes a predetermined matching threshold, where such 

70 threshold Is passed by such reference synthetic seismogram as providing the best match (i.e. highest 
crosscorrelation value) to the reference seismic trace, T28. 

Refenring again to FIG. 3, the seismic wavelet as derived above is indicated at 38. the convolution 
operation is represented by and the resulting reference synthetic seismogram is shown at 40. 
Reference synthetic seismogram 40 is also shown in FIG. 2 as being superimposed over reference T28. 

75 The following portions of the text entitled Seismic Stratigraphy , vol. 9, edited by Bob Hardage, 
Geophysical Press, London- Amsterdam, 1987, are referenced, and hereby incorporated by such reference, 
with respect to background information well known to those skilled in the art: seismic wavelet derivation 
techniques, pp. 257-268; production of synthetic seismograms, pp. 74-76; and wavelet shaping and 
processing, pp. 238-257. 

20 ' 

5. Deriving Set of Modified Log-Pairs 

A set of a predetermined number of modified log-pairs are. derived in accordance with this step. Each 
modified log-pair is different from one another and comprises a velocity log and a density log, each of 

25 which logs is a modified version of the respective velocity log and density log of the reference log-pair. 
Each of the modified log-pairs con-esponds to a petrophysical property or properties also associated with 
the reference log-pair, but the value(s) of such property or properties are different than the known value(s) 
associated with the reference log-pairs. The modified log-pairs should preferably be representative of 
probable and reasonable variations of the petrophysical property or properties for the various nonreference 

30 lateral locations of layer 12. 

According to a preferred embodiment and with reference to FIGS. 4A-4C, examples of possible logs in 
accordance with the Invention are shown as being a function of depth or time (i.e. two-way travel time 
obtained by conventional depth to time conversion), where such depth or time increases in the direction of 
the illustrated an-ows and where Vsh. Sw, v, and p are scaled according to typical values. Vsh, *, and Sw 

35 are in percent, v is in feet/second, and p is in g/cm^ . 

FIG. 4A illustrates Vsh, *, and Sw logs 26, 28, and 30, respectively, which are representative of 
petrophysical properties associated with the reference log-pair comprising velocity log 32 and density log 
34. Other useful petrophysical properties, also being associated with the reference log-pair, can be derived 
from Vsh log 26. $ log 28, and Sw log 30 and are described below. 

40 Pay (as previously defined) intervals, indicated in black in FIG. 4A, are those intervals of depth or time 
for which Vsh log 26, ^ log 28, and Sw log 30 exceed predetermined minimum "pay threshold" values. In 
terms of the depth dimension, several petrophysical properties can be defined as follows. The "gross pay 
thickness" is defined as the total thickness between the upper and lower limits of pay, which is less than or 
equal to the thickness of layer 12 at the reference lateral location, L28. The "net pay thickness" is the 

45 combined thickness of only the pay intervals. The "net pore thickness" is the sum of a series of * X AD 
products, where each such product corresponds to a different depth increment AD within a pay inten^al and 
between the upper and lower limits of pay, * is the porosity associated with depth increment AD, and 
where there are a predetermined number of depth increments between the upper and lower limits of pay 
which is generally equivalent to the depth increments at which log data was collected in step 2. The 

50 "hydrocarbon pore volume" is the sum of a series of * X (100%-Sw) X AD products, where ^ and AD are 
as defined with respect to net pore thickness and where Sw is the water saturation (in percent) associated 
with AD. The "net porosity" is the sum of a series of * values corresponding to respective depth 
increments AD, divided by the number of depth increments AD, where * and AD are as defined with 
respect to net pore thickness. In other words, net porosity is the average porosity in the pay intervals 

55 between the upper and lower limits of pay, and is equivalent to net pore thickness/net pay thickness. 

FIGS. 4B and 4C illustrate examples of two modified log-pairs, each comprising a velocity log and 
density log, and also examples of Vsh, ^> and Sw logs from which the modified log-pairs were derived in the 
manner discussed previously. It can be seen that the Vsh,*. and Sw logs in FIGS. 4B and 4C have been 
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modified with respect to amplitude as compared to the corresponding logs in FIG. 4A. Rather than 
modifying Vsh, *, and Sw concurrently as in FIGS. 4B and 4C, only one or a combination of only two of such 
properties could be modified in the derivation of the corresponding modified log-pairs. It should be apparent 
that layer thickness could also be modified. 

Each of the modified log-pairs resulting from this step have associated therewith the petrophysical 
property value(s) which were subject to modification in deriving the modified log-pairs. Such value(s) 
preferably correspond to a petrophysical property or properties which have a single recordable or storable 
value to be associated with any one modified log-pair, such as the pay thickness, net pore thickness, 
hydrocarbon pore volume, and/or net porosity which are derivable (as discussed above) from Vjh, and/or 
Sw togs. Or, such desired single values as associated with corresponding modified log-pairs could be 
average or mean Vsh, *. and/or S^ values derivable from the corresponding logs. 

The number of modified log-pairs provided in this step can vary widely, depending upon the known 
variability of the geology of the area being investigated and the areal extent of the layer of interest for which 
petrophysical property values are desired (I.e. area defined by the nonreference lateral locations). Typically, 
about 10 to about 200 modified log-pairs are derived according to the illustrated embodiment employing a 
single reference lateral location. 

As stated above, the modified log-pairs resulting from this step are preferably representative of 
probable variations of the petrophysical property value(s) over the areal extent of the layer of interest for 
which petrophysical property values are desired (i.e. area defined by the nonreference lateral locations). 
The extent to which modifications discussed above are permitted Is a somewhat subjective judgement by 
the geoscientist having knowledge of the geology of the area, so that constraints (minimums and/or 
maximums) upon modification actually made by a suitable computer program can be provided by the 
geoscientist as inputs to such program. 

6. Determine Modified Reflection Coefficients 

The velocity log and density log of each modified log-pair are employed to determine a corresponding 
modified reflection coefficient series (in the same manner as described in step 3 with respect to the 
reference reflection coefficient series), thereby resulting in a number of modified reflection coefficient series 
equivalent to the above-mentioned predetermined number of modified log-pairs. • 

7. Produce Modified Synthetic Selsmograms 

The seismic wavelet obtained In step 4 is convolved with each of the modified reflection coefficients (in 
the same manner as described in step 4 with respect to the reference synthetic selsmogram) to produce a 
con-esponding modified synthetic seismogram, thereby resulting in a number of modified synthetic seis- 
mograms equivalent to the above-mentioned predetermined number of modified log-pairs. Each of the 
modified synthetic selsmograms have reflection events corresponding to the upper and lower boundaries of 
the layer and also have associated therewith the value(s) of the petrophysical property or properties 
associated with the con^espondlng modified log-pair. 

8. Compare Nonreference Seismic Traces to Modified Synthetic Selsmograms 

A comparison window of each nonreference seismic trace Is compared to a corresponding comparison 
window of each of the modified synthetic seismograms, where each comparison window includes the 
reflection events corresponding to the upper and lower boundaries of layer 12. Each comparison window is 
preferably Identical in terms of time or depth and Is suffidently large to include the above-mentioned pair of 
reflection events for each trace. Preferred comparison techniques are discussed below, but any comparison 
technique can be employed according to certain broad aspects of the invention. 

Comparisons of this step can employ simple crosscorrelation so that each comparison yields a 
crosscorrelation value. Or, the comparisons can employ calculation of difference mismatch error (previously 
defined). 

Most preferably, however, a combination of crosscorrelation and difference mismatch error can be 
employed in the following series of steps as applied to comparison of a nonreference seismic trace and a 
modified synthetic seismogram, where each of the nonreference seismic trace and modified synthetic 
seismogram is defined by a series of data points at predetermined sample intervals (I.e. time): (i) shifting 
the modified synthetic seismogram with respect to the nonreference seismic trace a predetermined number 
of times with different corresponding sample shifts; (li) crosscorrelating the modified synthetic seismogram 
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and the nonreference seismic trace for each of the shifts in (i) to give a crosscorrelation value for each shift; 
(ill) selecting the shift which gives the maximum crosscorrelation value; and (iv) calculating the difference 
mismatch error between the modified synthetic seismogram, as shifted in accordance with the shift selected 
in (Hi), and the nonreference seismic trace. Of course, it is also possible in the above steps to shift the 
5 nonreference seismic trace or to shift both the nonreference seismic trace and modified synthetic 
seismogram. 

9. Select Modified Synthetic Seismograms for Each Nonreference Seismic Trace 

10 With respect to each nonreference seismic trace (T1-T27 and T28-T55), those modified synthetic 
seismograms are selected which match such nonreference seismic trace sufficiently In step 8 to pass at 
least one predetermined matching threshold. 

Where the comparison in step 8 is by either crosscorrelation or difference mismatch error, the matching 
threshold could be designated as a minimum crosscorrelation value (i.e. 0.90 or 0.95) in the case of 

75 crosscorelation or a maximum difference mismatch error value in the case of difference mismatch error. 
The matching threshold could also simply be designated as being passed by a preselected number n of 
modified synthetic seismograms whose comparison with a nonreference seismic trace provides the top n 
matches (i.e. n modified synthetic seismograms of the total number of such seismograms having the 
corresponding n highest crosscorrelation values or n lowest difference mismatch error values). 

20 Where the comparison in step 8 is by a combination of crosscon'elation and difference mismatch error, 
either or both of the above-described matching thresholds associated with crosscorrelation and difference 
mismatch error could be employed. The comparison procedure used in a subsequent example utilizes a 
matching threshold with respect to crosscorrelation as well as a matching threshold with respect to 
difference mismatch error. 

25 

10. Assigning Values of Petrophysical Property or Properties to Each Nonreference Seismic Trace and 
Corresponding Nonreference Location 

To each nonreference seismic trace there is assigned the values of the petrophysical property or 
30 properties associated with the modified synthetic seismograms which are selected in step 9 with respect to 
such nonreference seismic trace, thereby providing a range of such values associated with such non- 
reference seismic trace and with the corresponding nonreference lateral location. By way of example with 
respect to the petrophysical property net porosity and a particular nonreference lateral location, this step 
would result In a range of possible net porosity values corresponding to the nonreference lateral location, 
35 where such range consists of a highest possible net porosity value, a lowest possible net porosity value, 
and the most likely net porosity value which is associated with the modified synthetic seismogram best 
matching the nonreference seismic trace corresponding to the nonreference lateral location. 

The values assigned as discussed above can be displayed in any convenient manner, such as a 
numerical display of the lowest, highest, and most likely value(s) for each nonreference lateral location, a 
40 graphical plot display of the various values for each nonreference lateral location, a color coded map of 
most likely values, etc. 

Example 

45 This example demonstrates the effectiveness of the invention in determining a petrophysical property 

associated with a subterranean layer. 

Steps 1-10 of the inventive method as previously described were carried out with respect to a field of 

over 3,000 acres and a particular layer of interest In such field having upper and lower boundaries at depths 

of about 10,200 feet and 10,325 feet, respectively. Rfteen reference wells at 15 corresponding reference 
50 lateral locations in the field were employed to determine an average value of a petrophysical property (in 

this case, net porosity) associated with the layer for each of four nonreference wells at four corresponding 

nonreference lateral locations. Details of each step are given below. 

1. Reference seismic traces corresponding to the reference wells and nonreference seismic traces 
corresponding to the nonreference wells were taken from a set of 3-D seismic data for the field. Each of 

55 such seismic traces corresponded to the lateral locations of the wells or to the locations closely adjacent 
(i.e. within 25 feet) to the well location, 

2. Velocity and density logs for each of the reference wells were obtained as follows for the layer of 
interest, where each pair of such velocity and density derivative logs make up a reference log-pair. 
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Velocity (sonic) and density wireline well logs were employed in combination with gamma, spontaneous 
potential, and resistivity wireline well logs to derive a porosity (*) log. A water saturation (Sw) log was 
derived from resistivity or induction wireline wall logs. A lithological (Vsh) log was derived from a gamma 
wireline well log. The thus derived porosity, water saturation, and lithological logs were employed to 

5 derive velocity (sonic) and density logs via standard transit time (AT, travel time of a seismic wave 
through 1 foot of the layer, equivalent to the inverse of velocity) and density (p) equations, as are set 
forth as equations (1) (Wyllie time average equation) and (2) (bulk density equation) in the article 
"Incremental Pay Thickness Modeling of Hydrocarbon Reservoirs" by Dennis B, Neff, Geophysics , vol. 
55, no. 5 (May 1990), which is hereby incorporated by reference. 

70 3. A total of 15 reflection coefficient series were derived from the 15 corresponding reference log-pairs. 
4. Six different test seismic wavelets were extracted from a 500 millisecond window of a line of 50 traces 
taken from the above-mentioned set of 3-D seismic data, wherein each of the six test wavelets were 
based on different amplitude, frequency, and/or phase values. Each of such test wavelets were 
convolved with the reference reflection coefficient series of five selected reference wells to thereby 

75 produce five reference synthetic seismograms for each test wavelet. Each set of five reference synthetic 
seismograms was compared to the corresponding reference seismic traces. Such comparison was 
Interpretative, or subjective, In nature, relying primarily on crosscorrelation. One of the test wavelets was 
selected, based on the above-mentioned comparison, as the wavelet producing the reference synthetic 
seismograms best matching the corresponding reference seismic traces, 

20 5. For each of the 15 reference wells, between 25 and 50 modified log-pairs were derived by varying 
porosity and/or water saturation, as well as layer thickness, with respect to the known porosity and water 
saturation logs from step 2. Variations of amplitudes of the known porosity and/or water saturation logs 
ranged from about 10-25%, and variations in layer thickness were made In 20 foot increments of not 
more than ± 40 feet. These variations were made in a relatively random manner, varying thickness in 

25 each case, sometimes varying porosity and water saturation concun'ently, and in other cases varying 
either porosity or water saturation only. A total of 457 modified log-pairs resulted, where each such log- 
pair has associated therewith a lithological composition log. a porosity log, and a water saturation log 
from which net pay thickness and net pore thickness values are determined in the manner described in 
the detailed description of step 5 by assuming minimum "pay threshold" values. 

30 6. A modified reflection coefficient series was determined for each of the 457 modified log-pairs to 
thereby result in 457 modified reflection coefficient series. 

7. Each of the 457 modified reflection coefficient series was convolved with the seismic wavelet selected 
in step 4 to result in 457 modified synthetic seismograms, where each such modified synthetic 
seismogram has associated therewith a net pay thickness value and net pore thickness value from step 

35 5. 

8. A comparison window (40 milliseconds) of each nonreference seismic trace from step 1 was 
compared to a corresponding comparison window of each of the 457 modified synthetic seismograms by 
a combination of crosscorrelation and difference mismatch error. The upper and lower limits of each 
comparison window included reflection events corresponding to the respective upper and lower reflecting 

40 boundaries of the layer, such that such reflection events were centered within such comparison window. 
The comparison was carried out In accordance with the prefen-ed comparison procedure comprising 
substeps (i)-(iv) as described above in the detailed description of step 8, wherein 1 1 sample shifts were 
employed in substep (i). Consequently, the comparison of each nonreference seismic trace to each of 
the 457 modified synthetic seismograms resulted in a maximum crosscorrelation value and difference 

45 mismatch error value corresponding to each modified synthetic seismogram. 

9. For each nonreference seismic trace, 15 of the 457 modified synthetic seismograms were selected as 
having the 15 highest crosscorrelation values (passing a first matching threshold). Of the thus selected 
15 modified synthetic seismograms, seven were selected as having the seven lowest difference 
mismatch error values (passing a second matching threshold). Such seven modified synthetic seismog- 

50 rams are hereafter denoted as modified synthetic seismograms 1-7, where seismogram 1 has the lowest 
corresponding difference mismatch error value and seismogram 7 has the highest corresponding 
difference mismatch eror value. 

10. For each nonreference seismic trace, the net pore thickness and net pay thickness values associated 
with the selected modified synthetic seismograms 1-7 are assigned to such nonreference seismic trace, 

55 to thereby provide a range of net pore thickness and net pay thickness values associated with such 
nonreference seismic trace and its conrespondlng nonreference well. 

To facilitate a determination of the effectiveness of the invention, a single net porosity value (an 
average) was calculated for each nonreference well based on net porosity values obtained by the invention, 
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and a comparative net porosity value was obtained for each nonreference well by conventional means for 
the purpose of comparison. 

A single net porosity value (average) was determined for each nonreference well in accordance with the 
invention by first calculating a net porosity value for each of the corresponding modified synthetic 
seismograms 1-7 from the net pore thickness value and net pay thickness value (net pore thickness/net pay 
thickness) associated with each such seismogram, and then determining an average of the thus calculated 
net porosity values 1-7 (NPV1-NPV7) corresponding to modified synthetic seismograms 1-7. The average is 
calculated from the following formula: 

[(3 X NPV1) + 2X (NPV2 + NPV3 + NP4) + NPV5 + NPV6 + NPV7]/12. 

Comparative net porosities were obtained for each of the nonreference wells as follows, employing data 
obtained by wireline well logs of the layer of Interest. For each nonreference well, lithologlcal (Vsh), porosity 
(*), and water saturation (Sw) logs were derived in the same general manner as such logs were derived for 
the reference welis. From these logs, the pay inten^als were determined assuming the same "^pay 
threshold" values employed in step 5 of the invention as carried out in this example. A net porosity value 
was determined from the porosity log. assuming the previously determined pay intervals, by calculating the 
average porosity in the pay intervals in the manner discussed In the detailed description of step 5 (sum of 
porosity values at predetermined increments within pay intervals, divided by number of increments). 

The Table illustrates the above net porosity (*) results so as to provide a clear comparison between the 
net porosity values obtained by the invention and the comparative net porosity values. The Table sets forth 
the net porosity values, the variance of the Invention net porosity from the comparative net porosity, as well 
as the en'or (absolute value of the variance/comparative net porosity. 

TABLE 



Well 


Comparative Net $ (%) 


Invention Net * (%) 


Variance (%) 


Error (%) 


1 


23.1 


21 


-2.1 


9 


2 


29.0 


25 


-4.0 


14 


3 


15.3 


18 


+ 2.7 


18 


4 


25.3 


22 


-3.2 


13 



The Table clearly shows the excellent accuracy of the invention in determining a petrophysical property 
of a subterranean layer at a particular lateral location, based on limited log data from only 15 other lateral 
locations and also seismic data. It should be noted in particular that the Invention could be similarly 
employed to determine a petrophysical property or properties associated with the layer at any lateral 
location in the field of this example, which as noted above covers over 3,000 acres. In effect, the invention 
Integrates limited log data and seismic data for a layer of interest in a particular field so as to enable fast 
and economical determination of a petrophysical property or properties of such layer at any lateral location 
in the field. 

Computer Program 

Five important subroutines of a computer program for accomplishing data processing steps of the 
invention are set forth in Appendix I. Such subroutines are written in "C" language for a Sparc-10 computer 
manufactured by Sun, and is self explanatory to one skilled in the use of the Sparc-10 computer. 

"Subroutine A" generates a set of modified log-pairs (each comprising a 1/v (AT) log and a p log), and 
requires as input data a reference set of logs in digital form, including Vs^, *, and Sw logs. "Subroutine B" 
reads into computer memory a set of nonreference seismic traces, tfie estimated time position on such 
traces corresponding to the upper boundary of the layer of interest, and ttie modified synthetic seismog- 
rams resulting from convolution of a seismic wavelet with modified reflection coefficients determined from 
the modified log-pairs. "Subroutine C" calculates gross pay thickness, net pay thickness, net pore 
thickness, net porosity, and hydrocarbon pore volume values for each of the modified log-pairs and 
coHBsponding modified synthetic seismograms. and requires as Input data each modified set of Vsh, and 
Sw logs from which each of the modified log-pairs was derived. "Subroutine D" compares each non- 
reference seismic trace to each modified synthetic seismogram by crosscorrelation so as to calculate a 
crosscorrelation value for each of a number of shifts of a particular modified synthetic seismogram. and 
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requires a desired number and magnitude of slitfts as input data. "Subroutine E" determines the maximum 
crosscorrelation value for each modified synthetic seismogram as compared to a nonreference seismic 
trace, and calculates the difference mismatch error for each modified synthetic seismogram as shifted to 
give the maximum crosscorrelation value. "Subroutine F" selects m modified synthetic seismograms as 
compared to a nonreference trace having the m highest maximum crosscon'elation values, and of such m 
seismograms selects n seismograms having the n lowest difference mismatch error values, where m and n 
are integers and n < m. Subroutine F therefore requires m and n as input data. Subroutine F also provides 
an output of the gross pay thickness, net pay thickness, net pore thickness, hydrocarbon pore volume, and 
net porosity values con-esponding to each of the selected m seismograms, and calculates and provides as 
an output a weighted average of such values for the selected n seismograms. 

Conclusion 

Thus, there is provided by the present invention an effective method of determining the value of a 
desired petrophysical property of a layer at any lateral location thereof which requires a minimal amount of 
log data. Obviously many modifications and variations of the present invention are possible in light of the 
above teachings. It is therefore to be understood that within the scope of the appended claims the invention 
may be practiced othenwise than as specifically described. 
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SUBROUTINE A 



int compute num states (layrs) 
struct layr^Tnfo *layrs; 

{ 

int tot^states, i, j; 

/* Calculate number of steps per layer, number of states per layer, */ 
/* partial step sizes for each layer, and total number of states for */ 
/* this model. */ 

tot_states « 1; 

for(i«0; i<layrs->nl; i++) /* For each layer */ 

20 layrs->nsteps[i] » (int) (layr8->nsaap(i] / layrs->step_si2e[i] ) ; 

if ((layrs->nsampti] % layr8->step_size(i] ) 1-0) 

layrs->last_step size[i] » layrs->nsanp[i] - (layrs->nsteps[i] * 

layrs->step slze[i]); 
layrs->nsteps [T]++ ; 

) 

else 

layrs->last_step_si2e[i] = layrs->step_sizeCi] ; 

if ((layrs->inc_typeCi] 0) || (layrs->inc_type[i] i) ) 

layrs->nstates[i] » layrs->nsteps[i] + 1; 
if (layrs->inc_type[i] «- -1) 
( 

layrs->nstates [ i ] «» 2; 
for(j«l; j<layrs->nsteps[i) ; j++) 
layrs->nstates[i] *« 2; 

tot_states *= layrs->nstates[i] ; 

35 ) ~ 

/* Calculate layer divisors. */ 

layrs->d[0] = tot_states / layrs->nstates [0] ; 
for(i=l; i<layrs->nl; i++) 

layrs->d[i] = layrs->d[i-l] / layrs->nstates[i] ; 



return (tot_states) ; 



void compute_layer_8tates( layrs, mod^stata) 
struct layr_info *layrs; ~ 
int mod state; 



( 



55 



int j; 

for(j=0; j<layrs->nl; /* Update current states */ 

( /* of all layers ► */ 

layrs->stateC j) « (int) (mod_state / layrs->d( j] ) ; 

mod_Btate — (layrs->state[ j] * layrs->d[ j]) ; 

> 
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void bulld_sngl coap^log (layrs, init, repl, conp) 
struct layr Info *layrB; 
float *init7 *rcpl, *conip; 



{ 



int j, k, offsets- 
unsigned short bitpat; 

for(j=0; j<layrs->nl; /* ^or each layer */ 

^ switch ( layrs->inc_type [ j ) ) 

^ case 0: /* " top-down */ 

if (layrs->state[j] 0) , 

copy(init + layrs->l_offset[j], comp + layrs->l_of fset[ j] , 
layrs->nsamp[ j ] ) ? 

else 

^ if (layrs->stateCj] ~ layrs->nsteps[ j ] ) 

copy{repl + layrs->l_offset[ j ) , comp + layrs->l_offsetC j] , 
layrs->nsanpt j ]T; 

else 

copy(repl layrs->l_of f setC j ] , comp + layrs->l_off set[o ] , 
(layrs->step_si2e[j) * layrs->stateC j ] ) ) ; 

} 

break 7 

case 1: /* If bottom-up.... */ 

if {layrs->stateCj] 0) 

copy(init + layrs->l_of fsetC j ] , comp + layrs->l_of fsctC j] , 
1 ayrs->nsamp [ j ] ) ; 

else 

^ offset = layrs->l_offsetCj] + layrs->nsampC j ] - 
(layrs->state[j) * layrs->step_sizeC j ] ) ; 
if (layrs->state[j3 -= layrs->nsteps[ j ] ) 

copy(repl + layrs->l_offset[ j] , comp + layrs->l_offsetC j ] , 
layrs->nsampC j ] ) ; 

else 

copy (repl + offset, comp + offset, (layrs->step_si2e[ j ] * 
layrs->state[ j ] ) ) ; 

) 

break; 

case -1: /* If random */ 

bitpat = (unsigned int) layrs->state [ j ] ; 
for(k=0; k<layrs->nsteps[j] - 1; k++) 

^ offset « layrs->l_offset[j] + (k * layrs->step_size[ j ] ) ; 
. if (bitpat & 0x0001) 

copy (repl + offset, comp + offset, layrs->step_size[ j ) ) ; 
else 

copy(init + offset, comp + offset, layrs->step_sizeC j } ) ; 



bitpat »- 1; 

) 

offset +- layrs->step_8ize( j ] ; 
if (bitpat & 0x0001) 

copy(repl + offset, comp + offset, layrs->last_step_si2e[ j ] ) ; 
else 

copy(init + offset, comp + offset, layrs->last__step_size[ j]) ; 
break; ~ 

} 

) 

) 
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void build_cuinu coinp_log(layrs, init, repl, comp) 
struct layr_Tnfo *layr3; 
float ♦init, *repl, *comp; 



( 



int j, )c, offset; 

unsigned short bitpat, prev_bitpat, ord; 



for(j«o; j<layrs->nl; /* For each layer */ 

{ 

if (layrs->state[j] != layrs->prev_state[ j ] ) /* If state has changed.. ♦/ 
( 

switch (layrs->inc_type[j ] ) 
( 

case 0: /* If top-down.... */ 

if (layrs->state[j3 " 0) 

copy(init + layrs->l_of fsetC j] , comp + layrs->l_offset[ j ] , 
layrs->nsa»p £ j ] ) ; ~ 

else 

( 

offset « layrs->l_offset[j] 4 { (layrs->state[ j ] - 1) * 

layrs->step_si2eC j ] ) ; 
if (layrs->state[ j ] layrs->nsteps( j } ) 

copy (repl + offset, comp + offset, layrs->last_step_si2e[ j 
else 

20 copy (repl + offset, comp + offset, layrs->step_si2eC j ] ) ; 

) 

break; 

case 1: /* If botton-up . . • • */ 

if (layrs->stateCj] «« 0) 

copy(init + layrs->l_offset[j3 , comp + layrs->l_offsetC j ] , 
layrs->nsampc j ) J; 

else 
( 

offset = layrs->l_of fset [ j ] + layrs->nsampC j ] " 

(layrs->state[ j ] * layrs->step_size[ j] ) ; 
if {layrs->state[ j ] layrs->nsteps[ j ] ) 

copy (repl + layrs->l_of fset[ j ] , comp + layrs->l_offset[J ] , 
layrs->last_step_size[ j ] ) ; " 
else " " 

copy(repl + offset, comp + offset, layrs->step_si2eC j ] ) ; 

} 

break; 

case -1: /* If random */ 

bitpat = (unsigned int) layrs->state[ j ] ; 
35 prev__bitpat ■ (unsigned int) layrs->prev_state[j] ; 

ord = bitpat ^ prev_bitpat; 
for(k=0; k<layrs->nsteps[ j ] - 1; k++) 
( 

if (ord & 0x0001) 
( 

40 
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offset - layrs->l_offsetC j] + (k * layrs->step_si2e[ j ] ) ? 
if(bitpat & 0x0001) 

copy(repl + offset, comp + offset, layrs->step_si2e[ j ] ) ; 
eXse 

copy(init + offset, comp + offset, layrs->step_si2e[ j J ) ; 
bitpat »- 1; 
ord »- l; 

} 

> 

if{ord & 0X0001) 

^ offset » layrs->l_offsetCj3 + (k ♦ layrs->step_size( j ] ) ; 
if (bitpat & 0x0001) 

copy(repl + offset, comp + offset, layrs->last_step_size[ j j^* 

copy(inlt + offset, comp + offset, layrs->last step_si2e[ j 1 • 

} 

break; 
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SUBROUTINE B 

I include <atdio.h> 
♦include <lo.h> 
5 I include <nalloc.h> 
I include <Btring.h> 
I include <stdllb.h> 
# include "nui.h" 

void IbiaTorcee( float *, float *, int) ; 
int flipit(int) ; 
'0 long longflip(long) ; 

int get_num_traces(FILE *segfile, long nuabytes) 
( 

int numsamp =» 0, nuntr - 0; 
long loc » 3714; 

75 

while (loc < numbytes) 
( 

fseek(segf ile, loc, SEEK^SET) ; 
fread(&numsamp, sizeof(int), 1, segfile) ; 
numsamp « flip it (numsamp) ; 
nujiitr++; 

20 loc += (240 + (4 * numsamp)); 

} 

return (numtr) ; 

) 



25 



long get_si2e(FILE *pcfile) 
{ 

int fh; 

fh « fiXeno(pcfile) ; 
retum(filelength(fh) ) ; 

) 

30 void hpprint( struct modattrib *iptattrib, int first, int last) 
( 

FILE *prnfile; 
int i; 



35 



45 
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pmfile « fopen("prn", "w"); 

fprintf (prnfile, "\xlB&]c2S") ; 
fprintf (pmfile, "\)clB&19D") ; 



fprintf (prnfile, "IPT Attribute Printout\n\n\n") ; 

fprintf (prnfile, "CDP PK, AMP TRO AMP APP ISO ACT ISO TIM TRO TIM PK. ") ; 

fprintf (prnfile, "TIM TOP TIM BAS NET PAY NET PGR GRO PAY POR VOL »') ; 

^ fprintf (prnfile, " DT SAG TOP-PK BAS-PK BAS-TR\n") ; 

fprintf (pmfile, »• ") ; 
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fprintf(prnl!ile, " 'Z'JIZ vn"l • 

fprintf (prnfile, " 

for(i«first; i<"last 



( 



) 



fprintf (prnfile, 
fprintf (prnfile, 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
fprintf (prnfile, 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 
iptattrib 



%3d iptattribCi} .cdp_nun) ; 

%7.4f %7.4f %7,2f %7.2f *7.2f %7.2f %7»2f «, 

i] .pea)c_aiiip, 

i ) . trough amp , . 

i] .appar_Tso, 

i) .act_iso, 

i] .tine_trough, 

i] ,tinej?^^^f 

il .tittG top_pay) ; 

*7.2f %7,2f %7.2f %7.2f %7.2f *7.2f %7.2f *7.2f %7.2f\n« 
i) .time^base jpay, 
ij.netjay, 
ij .net_por_ft, 

i] ,gros3_j>ay, 
i] .hc_por_vol, 
i].dt_sag, 
i) .top_pea)c, 
ij .base_peaX, 
i) .base_trough) ; 



fprintf (prnfile, "NxlB") ? 
fprintf (prnfile, "E") ; 
fclose (prnfile) ; 



void get seg_tr(FILE *segfile, long stbyte, int nun^pts, float *output) 
C 

float *tempibm; 

tempibm » calloc(nuD_pt3, sizeof (float) ) ; 
fsee)c(segfile, stbyte, SEEK_SET) ; 
f read (temp ibn, sizeof (float) , num_j>ts, segf ile) ; 
IbmToIeee (tempibm, output, nunjts) ; 
free(teT&piba) ; 

) 

void get iptsegtrc hdr(FILE *segfile, long fileloc, struct nodattrib *output) 

( " 
union 

{ 

float fval(60]; 
long lvai(60]; 
int ival[120]; 
) trc_hdr; 
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f seek (segf lie, fileloc, SEEK_SET) ; 
fread(trc_Mr.fval, slzeof (float) , 60, segf ile) ; 



output - > s t a r t_by t e 
output- >cdp_nuia 
output- >nsainp 
output->delay 
output->firate 



- fileloc + 240; 

- longflip(trc_hdr.lval(5]) ; 
« fliplt(trc hdr.ivalC57]) ; 

- flipit(trc2hdr.ival[54]) ; 
= (fllpit(trc hdr.lval(58]))/1000; 



ouT:pu^-;>Bi:a\.e - \ i.j.At^-i.v v .* j ---- , 

IbmToIeeeC (float *) &trc_hdr.fval[43] , (float ♦) & output ->trough_ainp, 17); 
} 

void get_segtrc_hdr(FILE *segfile, long fileloc, struct trcattrib *output) 

union 
< 

float fvaiceo]; 
long lval(60]; 
int ival[120]; 
) trc_hdr; 



fseek(segfile, fileloc, SEEK^SET) ; 
fread(trc_hdr.fval, sizeof (float) , 60, segfile) ; 



) 



output->start_byte 
output->cdp_nujn 
output->nsaiap 
output ->de lay 
output->srate 
IbiaToIeeeC (float *) 



• fileloc + 240; 

- longflip(trc hdr. IvalCS] ) ; 
= flipit(trc hdr.ivalCS?]) ; 

- flipit(trc'"hdr,ivalC54]) ; 

» (flipit(trc hdr.ivalC58]))/10O0; 

&trc_hdr.fvaT[50], (float*) &output->top_pay_pick, 1) ; 



void get_l]a3csegtrc_hdr(FILE *segfile, long fileloc, struct trcattrib *output) 



{ 



union 

{ 

float fval[60} ; 
long lval[6oj; 
int ivalC120]; 
} trc_hdr; 

fseek(segfile, fileloc, SEEK_SET) ; 
fread(trc_hdr. fval, sizeof (float) , 60, segfile); 



=• fileloc + 240; 
« longflip(trc__hdr.lvalC5]) ; 
» flipit(trc3dr.ival[57]) ; 
» -flipit(trc_hdr.ival[52]) ; 
= (flipit (trc^hdr . ival [58] ) )/1000; 
IbniToIeee( (float *) &trc hdr.fval[50] , (float *) &output->top_pay_pick, 1) ; 
if ( (output->top_payjpick""> 6000) || (output->top_pay_pick <« 0)) 



ou t put-> St a r t_by t e 
output->cdp_nuia 
output->nsamp 
output->delay 
output ->srate 
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output->error_flag - 1; 
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SUBROUTINE C 



void payatats (model, layers, trstats, paylog, porlog, swlog) 
IPTModelKeader *]nodel; 
IPTLayerlnfo *layer8; 
IPTTraceParaaa ♦trstats; 
float *paylog; 
float *porlog; 
float *swlog; 

{ 

float *payptrl, *payptr2, *porptrl, *porptr2, *swptrl, *swptr2; 
int i; 

for(i-0; i<iftodel->nuin layers; 
{ 

payptrl = paylog + layers[i} .layer_of fset; 
payptr2 = payptrl + layers [i] .nsamps^layer - 1; 
porptrl « porlog + layers [ i] .layer_off set; 
porptr2 = porptrl + layers [i] .nsamps^layer - 1; 
swptrl = swlog + layers [ i] .layer_off set; 
svptr2 = swptrl + layers [ i] .nsamps^layer - 1; 

trstats->net_pay[i] = 0; 
trstats->layer_pore_feet[i] = 0; 
trstats->gross_pore^feet[l] » 0; 
trstats->net_pore_feet[i] « 0; 
trstats->layer_hpv[i] « 0; 
trstats->gross_hpv[i] « 0; 
trstats->net_hpv[i] » 0; 

trstats->thicJcness(i] « payptr2 - payptrl + 1; 

while ( (*payptrl ~ 0) && (payptrl <= payptr2)) 
( 

trstats->layer_j>ore feet[i] *porptrl / 100; 

trstats->layer_hpv[T] += *porptrl * (100 - *swptrl) / 10000; 

payptrl++; 

porptrl++; 

swptrl-i-+; 

> 

while ( (*payptr2 == 0) && (payptr2 >= payptrl)) 
( 

trstats->layerjpore feet(i] += *porptr2 / 100; 
trstats->layer_hpv[T] *porptr2 * (100 - *swptr2) / 10000; 
payptr2 — ; 
porptr2 — ; 
swptr2 — ; 

) 

trstats->gross_pay(i] = payptr2 - payptrl + 1; 

if (trstats->grossjpay[i] » 0) 

( 

layers ( i ] . top_gross_of f set « { int) NULL_VAL; 
layers [ i] .base_gross_off set ■ (int) NULL^VAL; 
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) 

else 

^ layers [ i] .top_gross_off set = payptrl - paylog; 
layers [ i] .base_gross_off set = payptr2 - paylog; 

for(; payptrl<=payptr2 ; payptrl++, porptrl++, swptrl++) 

^ trstats->layer_pore feet[i] += *porptrl / 100; 

trstats->layer_hpv[T] *porptrl * (100 - *swptrl) / 10000; 
trstats->gross_pore feet[i] *porptrl /lOO; 
trstats->gross_hpv[T] +- *porptrl * (100 - *swptrl) / 10000; 
if(*payptrl l» 0) 

tr stats->net_pay [ i ] ++ ; 

trstats->netjore feet[i] += *porptrl / 100; 
trstats->net_hpv[I] += *porptrl * (100 - *swptrl) / 10000; 

) 

if (trstats->gross_j)ay[i] » 0) 
{ 

trstats->net_gross_ratio[i] « 0; 
trstats->gross_avg_porosity[i] « NULL^VAL; 
trstats->net_avg_porosity[i] = NULL^VAL; 
trstats->gross_avg_sw[i] » NULL_VAL; 
trstats->net_avg_sw[ij « NULL_VAL; 

) 

else 

trstats->net_gross_ratioCi] =* (float) trstats->net_payCi] / 

(float) trstats->gross_pay[i] ; 
trstats->gross_avg_porosity[i] = trstats->gross_pore_feet[i] / 

trstats->gross_pay[ i] * 100; 
trstats->net_avg_porosity[i] » tr stats ->net_pore_feet[i] / 

trstats->net_pay[i] * 100; 
trstats->gross_avg_sw[i] » ( (trstats->gross__pore_feet [i] - 

trstats->gross_hpv[i]) / trstats->gross_pore_f eet[i] ) 

* 100; 

trstats->net_avg_sw[i] « { {trstats->net_pore_f eet [i] - 

trstats->net_hpv[i]) / trstats->net_pore feet[i]) * 100; 

} 

trstats->layer__avg_porosity[i] « trstats->layer_pore_feet[i] / 

trstats->thic)cness[i] * 100; 
trstats->layer_avg_sw[i] = ( {trstats->layer_pore_f eet [ i] - 

trstats->layer hpv[i]) / trstats->layer_pore_f eet [i] ) * 100; 
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SOBROUTIWE D 

i include <stdlo.h> 
linclude <io.h> 
I Include <inalloc.h> 
linclude <string.h> 
i include <aath.h> 
linclude <stdlib.h> 
linclude "nui.h" 
linclude <conio.h> 

Idefine PAUSE printf ("\nPress any key to continue \n") ; getch() 

void IbinToIeee( float *, float *, int) ; 

int flipit(int) ; 

long longf lip (long) ; 

int get_nu2n_t races (FILE *, long); 

long get size (FILE *) ; 

void hpprint( struct nodattrib *, int, int) ; 

void get seg tr(FILE *, long, int, float *) ; 

void BhortXcTfloat *, float *, int, int, float *) ; 

void get_iptsegtrc_hdr(FILE *, long, struct modattrib *) ; 

void get_lj!i)csegtrc_hdr(FILE *, long, struct trcattrib *) ; 

void get segtrc_hdr (FILE *, long, struct trcattrib *) ; 

void aax~and_lag( float *, int, float *, int *) ; 

double sunisq( float *, int); 

float difference (float float int); 

float dual_sign_bit( float *, float *, int); 

void correlate (float *, int, float *, int, float *) ; 

void hilbt oper( float *, int); 

void phs_and_lag( float *, int, float int, float *, float *, int *) ; 

void nuifunc(void) 
{ 

struct nodattrib *iptattrib; 

struct trcattrib *seisattrib; 

struct compare *conp; 

struct conpare^f ile_hdr *coiaphdr; 

char in_naitte(30], answer [10], trstrCS]; 

FILE ♦modfile, *scisfile, *compfile; 

int i, numtrm, numtrt, trace, index, fstpr, Istpr, numpts, stseunp, zerlag 

int delay, j, stpoint, k, diffpoint, startindex, endindex; 

int inincdpfil, maxcdpfil, mincdpproc, naxcdpproc, pad « 0, nlags; 

int fix_win_f lag, fix_win_len, skew_offset; 

long numbytesm, numbytest, fileloc, stbyte; 

float *inodseg, *trcseg, * crosscut, *hbtoper, *guadrature, *fullcross; 
float qfact, rel_seis_gain; 

/* Get name of SEGY format seismic trace file. V 

printf{"SEGY seismic trace file name? ") ; 
in_name [ 0 ] « • \o • j 
gets(in_name) ; 
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whllo((sol8file - fopen(ln_naino, "rb")) " NULL) 

^ printf ("\nFilG not foundl Try again. \n"); 
prlntf("SEGY seismic trace file name? ") ; 
5 in_name[0) - '\0' s 

gets(in_nane) ; 

} 

printf ("Checking file \n\n") ; 
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nuxnbytest « get sizeCseief ile) ; 

numtrt « get_n\m_traces(Beisf ile, numbytest) ; 

printf ("Number of bytes in trace file is %ld\n", numbytest)? 
printf ("Total number of traces in this file is *d\n\n", numtrt); 

printf ("Building trace attribute structure \n\n") ; 



seisattrib = calloc (numtrt, sizeof (struct trcattrib) ) / 
fileloc » 3600; 
if (IMKSOURCE) 

fGr(i»0; i<numtrt; i-H-) 

^ get lmJcsegtrc_hdr(seisfile, fileloc, &seisattrib(i] ) ; 

fileloc +« 240 + (4 * seisattrib ( i] .nsamp) ; 
} 

else 

for(i=0; i<numtrt; i++) 

^ get_segtrc hdr (seisf ile, fileloc, &seisattrib[i] ) ; 
25 fileloc +»"240 + (4 * seisattrib [ i] .nsamp) ; 

} 

/* Determine the sesimic CDP range to be processed. */ 



mincdpfil = (int) seisattrib [ 0] .cdp_nuB; 
maxcdpfil » mincdpfil + numtrt - 1; 



mincdpproc » mincdpfil; 

printf {"\nThe CDP range represented by your seismic file is %d to %d.\n\n", 

mincdpfil, maxcdpfil); 
printf ("Enter the first CDP to process 
printf ("Default = %dAn", mincdpfil); 
^5 printf («==«> "); 

answer(0] = '\0»; 

gets (answer) ; 

if (strlen (answer) 1= 0) 

{ 

mincdpproc = atoi (answer) ; 
40 while ((mincdpproc < mincdpfil) |1 (mincdpproc > maxcdpfil)) 
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) 



nincdpproc - mincdpfll; 

printf ("\nEntry out of range try agalnl\n"); 

printf ("~-> ") ; 
answerCO] - *\0' ; 
gets (answer) ; 
if (strlen (answer) J» 0) 

aincdpproc - atoi (answer) ; 

) 



maxcdpproc = maxcdpfil; 

printf ("\nThe CDP range represented by your seisaic file is %d to *d.\n\n", 

aincdpfil, aaxcdpfil) / 
printf ("Enter the last CDP to process. \n" ) ? 
printf ("Default « %dAn", naxcdpfil) ; 
printf ("-==> ") ; 
answerCO] « 'XO'; 
gets (answer) ; 
if (strlen (answer) 1« 0) 
( 

maxcdpproc - atoi (answer) ; 

while ((maxcdpproc < mincdpfil) || (maxcdpproc > maxcdpfil)) 
maxcdpproc » maxcdpfil; 

printf ("\nEntry out of range try again l\n"); 

printf ("===> "); 
answer! 0) = »\0' ; 
gets (answer) ; 
if (strlen (answer) i« 0) 

maxcdpproc » atoi ( answer) ; 

) 

} 

/* Determine the gain which should be applied to the seismic data to */ 
/* amplitude-balance it with the model data. */ 

rel_seis_gain = 1; 

printf ("\nEnter the relative gain factor to be applied to the seismic\n") ; 

printf ("data to amplitude-balance it with the model data.\n"); 

printf ("Default « %f.\n", rel_seis_gain) ; 

printf ("«==> "); 

answer (0) = '\0*; 

gets (answer) ; 

if (strlen (answer) 1= 0) 

rel_seis_gain = (float) atof (answer) ; 

while (rel_seis gain « 0) 

( 

rel_seis_gain « 1; 
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printf ("XnEntry out of range. ... .try againl\n"); 
prlntf("— > 
answer [0] - •\0'; 
gets (answer) ; 
if (strlen (answer) 1- 0) 

rel seis gain » (float) atof (answer) ; 

) 

> 

/* Get name of SEGY format IPT model traces file. */ 

prlntf ("\nIPT SEGY model file name? ") ; 
gets (in_name) ; 

while ({modfile « f open ( in_name , "rb")) " NULL) 

^ printf ("\nFile not foundl Try again. \n") ; 
printf("IPT SEGY model file name? "); 
in_name[0] » 'NO'; 
gets(in_name) ; 

} 

printf ("Checking file \n\n"); 

numbytesm = get^sizc (modfile) ; 

numtrm = get_niim_traces (modfile, numbytesm); 

printf ("Number of bytes in model file is %ld\n», numbytesm); 
printf ("Total number of traces in this file is %d\n\n", numtrm) 

printf ("Building model attribute structure \n\n") ; 

iptattrib « calloc (numtrm, sizeof (struct modattrib) ) ; 
fileloc = 3600; 

comphdr = calloc(l, sizeof (struct compare_f ile_hdr) ) ; 
comphdr->min_nt_pay = 5000; 
comphdr-'>min_nt_j)ore_ft = 5000; 
comphdr->min_hc_vol » 5000; 
comphdr •>max_nt3ay = 0; 
comphdr->max_nt_pore_ft = 0; 
comphdr->max_hc_vol « 0; 

for(i-0; i<numtrm; i++) 
{ 

get iptsegtrc hdr(modfile, fileloc, &iptattrib[i]) ; 
if (Iptattrib [T] » net jay < comphdr->min_nt_pay) 

comphdr- >min_nt_pay » Iptattrib [ i ] . net_pay ; 
if (iptattrib [i] .net_pay > comphdr- >max_nt_pay) 

comphdr- >max_nt_j)ay = iptattrib [ i] .net_pay; 
if (iptattrib [i] .net_por_ft < comphdr->min_nt_pore_ft) 

comphdr->min_nt_pore_ft = iptattrib [ i] .net_por_ft; 
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If (iptattribfl) .not_por_ft > comphdr->max_nt_pore_ft) 
comphdr->Bax ntj)ore_ft - iptattrib [ 1 ] . net jor_f t ; 

if (iptattrib[l)Thc_por vol < conphdr->»ln_hc_vol) 
conphdr->nin_hc vol""- iptattr ib [ i ) . hc^or^vol ? 

if (lptattrib[i).hcjpor_vol > comphdr->max_hc_vol ) 
coinphdr->nax_hc vol - iptattrib[i] .hc^or_.vol; 

flleloc +» 240 + (4 * iptattribCiJ.nsamp) ; 



/* Determine if fixed or variable windowing will be used. */ 



printf ("\nl8 the trace comparison window to be variable or fixed in lcngth?\ 
printf("V/F, Default « V\n") ; 
printf ("«=««> "); 
answer [0] = »\0'; 
^5 gets (answer) ; 

if (strlen (answer) =« 0) 
strcpy (answer , " v" ) ; 
if (strcmpi( answer, "f") >« 0) 
fix_win_flag * 1; 

/* If the window is to be variable, determine the pad and */ 
20 /* possible skew offset. */ 

if (Ifix win flag) 
( 

pad -"15; 

printf ("\nPlease enter the appropriate pad in &sec.\n"}; 
printf ("Default = 15 osec.\n") ; 
^ printf ("«==> ") ; 

answer [0] = 'XO' ; 
gets (answer) ; 
if (strlen (answer) !» 0) 
pad atoi (answer) ; 

30 skew^offset » 0; 

printf ("\nNonnally, the comparison window is centered on the model"); 
printf (" zone of interest. \n") ? 

printf ("Please enter any desired shift from this position in msec.\n"); 

printf ("Negative shifts window up, positive down.\n"); 

printf ("Default « 0\n") ; 

printf ("===> 

answerCO] « '\0'; 

gets (answer) ; 

if (strlen (answer) J« 0) 

skew^of f set = atoi (answer) ; 



40 /* If the window is to be fixed, determine the window length including V 
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/* pad and the offset between the model top and the window start. */ 

if (fix wln^flag) 
{ 

fix win_len - 40; 

prlntf ("\nPlease enter the length of the comparison window in nsec.\n") 

printf ("Default - 40 msecAn"); 

printf (••«=-> 

answer[0] - 'XO'; 

gets (answer) ; 

if (strlen( answer) I" 0) 

fix_win_len » atoi (answer) ; 

skew_offset « -15; 

printf ("\nPlease enter the the window start position in msec relative") 

printf (" to the top of the\nmodel zone of interest. \n") ? 

printf ("Negative means window starts above, positive below. \n"); 

printf ("Default « -15\n") ; 

printf ("«==«> ") ; 

answer [ 0 ] « • \0 • ; 

gets (answer) ; 

if (strl en (answer) 1» 0) 

s)cew_offset " atoi (answer) ; 



nlags = 11; 

printf ("\nFlease enter the number of lags at which to calculate a\n") ; 

printf (''correlation coefficient. Xn**) ; 

printf ("Default « ll\n") ; 

printf ("==> ") ; 

answer[0] = '\0'; 

gets (answer) ; 

if (strlen (answer) 1= 0) 

{ 

nlags » atoi (answer) ; 

if ((nlags % 2) = 0) /* Make sure nlags is odd. */ 

nlags++; 

) 

printf ("\nName for output statistics file? ") ; 
in^name [ 0 ] =» • \0 • ; 
gets(in_name) ; 

compflle = f open ( in_name , "wb") ; 

/* Fill remaining fields and write out compare file header. */ 

comphdr->st_cdp » mincdpproc; 
comphdr->cdp_inc = 1; 

comphdr->num_seis_trcs = maxcdpproc - mincdpproc + 1; 
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coaphdr->nun_modl_trc9 - nuatna; 
comphdr->xcorr_nlags - nlags; 
coittphdr->xcorr pad - pa^; 
comphdr->fix_wIn^flag - fix win flag? 
cQinphdr->fix win^len - flx^win len7 
conphdr->s)cew offset «> skew^offset; 
coaphdr->rol_sei9_galn - rel_seis_gain; 

fwrlte{comphdr, sizeof (struct coiapare_f ile^hdr) , 1, compf ile) ; 
free(conphdr) ; 

/* Display model attributes if desired. */ 

printf("\nDo you wish to view the IPT model attributes?\n«) ; 

printf("Y/N, Default « M\n»') ; 

printf ")7 

answer [0] - »\0«; 

gets (answer) ; 

if (strlen( answer) ~ 0) 

Btrcpy( answer, "no"); « 
if((Btrcmpi (answer, "y") — 0) 1 1 (strcmpi (answer, "yes") — 0)) 

* printf ("\nView attributes for which model trace (sequential from 1)? ") 
get8(trstr) ; 
trace - atoi(trstr) ; 
index - — trace; 



while (index>«0 && index<=numtrm-l) 

^ printf ("Starting byte 
printf ("COP number 
printf ("Number samples 
printf ("Delay 
printf ("Sample rate 
printf ("Peak amplitude 
printf ("Trough amplitude - %f\n", 
printf ("Apparent isochron - %f\n", 
printf ("Actual isochron « *f\n"- 



» %ld\n", 
» %ld\n", 
« %d\n", 
» %d\n", 
« %d\n\n" 
- %f\n» 



printf ("Trough time 
printf ("Peak time 
printf ("Time top pay 
printf ("Time base pay 
printf ("Net pay 
printf ("Net pore feet 
printf ("Gross pay 
printf ("Hyd pore volume 
printf ("DT sag 
printf ("Top peak 
printf ("Base peak 
printf ("Base trough 



- %f\n", 

- %f\n", 
= *f\n", 
» %f\n", 
» %f\n", 
« %f\n", 
= %f\n", 
« %f\n", 

- %f\n", 

- %f\n", 
= %f\n", 
« %f\n\n" 



iptattrib [ index 
iptattrib [ index 
iptattrib [ index 
iptattrib [index 
iptattrib ( index 
iptattrib [index 
iptattrib [ index 
iptattrib [index 
iptattrib [index 
iptattrib [ index 
iptattrib [ index 
iptattrib [ index 
iptattrib [ index 
iptattrib [index 
iptattrib [ index 
iptattrib (index 
iptattrib [index 
iptattrib [index 
iptattrib [index 
iptattrib [index 
iptattrib [index 



, start_byte) ; 
. cdp_num) ; 
.nsamp) ; 
•delay) ; 
. srate) ; 
,peak_amp); 
.trough amp) ; 
.appar_Tso) ; 
. act_iso) ; 
. time_trough) ; 
. time 3eak) ; 
. time_top_pay) ; 
. tioe_base_pay) 
. net_pay) ; 
. net_por_f t) ; 
. gross^pay) ; 
.hcjor_vol) ; 
. dt_sagT; 
. top jeak) ; 
,base_jpeak) ; 
. base_trough) ; 
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prlntf("View attributes for which model trace (sequential fron 1)? 

gets(trstr) ; 

trace • atoi(trstr) ; 

index - — trace; 

) 

) 

/* Print model attributes if desired. */ 

printf("\nDo you wish to print the IPT model attributes?\n") ; 
printf("y/N, Default - N\n»»); 
printf("=~> 
answer (0] » '\0*; 
gets (answer) ; 
if (6trlen( answer) ««=» 0) 
strcpy (answer, "no"); 
if ( (strcmpi (answer, "y") »» o) | | (strcmpi (answer, "yes") «= 0)) 

printf ("\nFirst model trace (sequential from 1) to print? ") ; 
scanf("%d", fifstpr) ; 

printf ("\n\nLast model trace (sequential from 1) to print? ") ; 
scanf("%d", filstpr) ; 
hpprint(iptattrib, fstpr, Istpr) ; 

) 

/* Print out processing information. */ 

if (fix_win flag) 

{ 

printf ("XnS elected windowing type is FIXED. \n**); 
printf ("The window length is %d msec.\n", f ix_win_len) ; 
printf ( "The window begins %d msec above the mSdel^pay 2one.\n\n", 
-s)cew_of fset) ; 

) 

else 

{ 

printf ("\nselected windowing type is VARIABLE. \n*') ; 

printf ("The window length includes a pad of %d msec above and below", 
pad) ; 

printf (" the model pay zone,\i^")' 

printf ("The window begins %d msec above the model pay zone.\n\n", pad 
skew_of f set) ; 

) 

hbtoper » calloc(HILBLN, sizeof (float) ) ; 
hilbt_oper(hbtoper, HILBIil) ; */ 

/* Loop through seismic traces, getting whole trace one<-at->a-time. */ 
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startlndex - mlncdpproc - (int) eeisattribCO] .cdp_nuia; 
endindex - startlndex maxcdpproc - aincdpproc + 1; 
comp - calloc(nuintrm, slzeof (struct compare)}; 

for (i«B tart Index; i<endindex; /* For each seismic trace */ 

^ if (Beisattrib[i] .error_flag «- 1) 

^ printf ("Missing horizon pick on seismic trace *ld. Trace s)cipped.\n"/ 
seisattribCi] .cdp^nua) ; 
for(j»0; j<numtnB; j++) 

^ coiap[j] .seis cdp « seisattrib[i] .cdp_num; 
comp[j] .mod_cdp « lptattrib[ j J .cdp^nun; 
comp C j ] . z er_xcorr • 0 ; 
comp C j ] . max_xcorr - 0 ; 
comp [j] .difference » 0; 
comp C j 5 . ph3_max_env - 0 ? 
comp ( j ] . quality"- 0 ; 
comp [ 5 ] . signbit " 0 ; 
comp [ j ] . lag_max_env « 0 ; 
comp [ j ] . meoc^lag « 0 ; 

fwrite(comp, sizeof (struct compare), numtrm, compf ile) ; 
continue; /* S)cip this trace if no horizon pick. V 

) 

printf ("Processing seismic trace %ld \n", seisattribCi] .cdp_nun) ; 

trcseg « calloc(seisattrib[i} .nsanp, sizeof (float) ) ; 

get_seg_tr(seisfile, seisattrib[i] .start_byte, seisattribCi) insanp, trcse 

Scale seismic up by user-supplied factor to match model amplitudes, */ 

for(j=0; j<seisattribCi] .nsamp; 
trcseg[j] rel_seis_gain; 

if(fix_wln flag) 

stpoint"- (int) (seisattribCi) .top^ay_pic)c + 0.5) 

- seisattribCi). delay + skew^offset - ( (nlags-l)/2) ; 

else 

stpoint « (int) (seisattribCi) .top_pay_plc)c +0.5) 

- seisattribCi). delay - pad + skew_offset - ( (nlags-l)/2) ; 

/* Loop through model traces, getting only desired portion. */ 

for(j=0; j<nunitrm; j++) 
{ 

delay = iptattribCj) .delay? 
if (fix win_flag) 
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numptB - fix win len + 1; 

stsamp - (int) {Tptattrib( j] . tiine_top_pay + 0.5) - delay 
skew offset; 

> 

else 

^ numpts - ((int) (iptattribf j ] . tiine_base_pay +0,5) - 

(int) (iptattribt j] .tiae_top_pay + 0.5)) + 2*pad + 1; 
stsamp - (int) (iptattrib[ j ] .tine_top_pay + 0.5) - delay - 
pad + skew^of fset; 

modseg *> calIoc(nuiapts, sizeof (float) ) ; 

stbyte - iptattrib[j).start_byte + (stsamp * 4); 

get_seg_tr(ttodfile, stbyte, numpts, modseg); 

/* Generate array of normalized correlation coefficients */ 

/* corresponding to a few lags above and below aligned nodel */ 

/* and seismic tops. */ 

crosscut B calloc(nlags, sizeof (float) ) ; 

shortXC(&trcseg[stpoint] , modseg, numpts, nlags, crossout) ; 

Load results of this comparison into compare structure. */ 

comp[ j ] .seis_cdp « seisattrib[i] .cdp_num; 
comp [ j ] • mod_cdp = iptattrib ( j ] . cdp_nua ; 
compCj] .zer_xcorr - crossout ( (nlags-l)/2] ; 

nax_and_lag( crossout, nlags, &comp[j3 .max^xcorr, 6comp[ j ) .max_lag) ; 

/* Compute full cross-correlation of windowed trace and model */ 
/* segments which resulted in best correlation coefficient above. */ 

fullcross » calloc( (2*numpts)-l, sizeof (float) ) ; 

correlate (&trcseg[stpoint + (nlags-l)/2 + comp[ j ] .max^lag) , 

numpts, modseg, numpts, fullcross) ; ^ */ 

/* Determine relative lag and instantaneous phase at */ 
/* peak of envelope of full cross-correlation. */ 

quadrature » calloc ( (2*numpts)+HIIiBLN-2 , sizeof (float) } ; 
phs_and_lag( fullcross, (2*numpts)-l, hbtoper, HILBLN, quadrature, 
&comp[ j] .phs_max_env, &comp[j] .lag_max_env) ; */ 

/* Calculate "quality factor" consisting of correlation */ 

/* coefficient downgraded by a factor based on the phase of */ 

/* the cross-corr, an indicator of its "skewness**. */ 

qfact = (PHSBASE - (ABS (comp[j] .phs_inax_env) ) ) / PHSBASE; 
if {qfact<0) 
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compCj] .quality - compC j] .aax_xcorr * qfact; */ 

free (crosscut) ; 
/* free ( f ullcross) ; 

free ( quadrature ) ; */ 

diffpoint - stpoint + coinp[ j ] .nax^lag + ( (nlags-l)/2) ; 

compel ] .difference - difference (fitrcseg [diffpoint ] , aodseg, numpts) 

conp[j).signbit « dual_8ign_bit (fctrcseg [diffpoint] , modseg, numpts) 

free (nodseg) ; 

fwrite(comp, sizeof (struct compare), nuntrm, compf ile) ; 
free(trcseg) ; 

) 

frce(hbtoper) ; */ 
f close (modf ile) ; 
fclose(s6isfile) ; 
fclose( compf ile) ; 
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SUBROUTINE E 



linclude <stdio.h> 
linclude <io.h> 
linclude <malloc.h> 
linclude <string.h> 
linclude <stdlib.h> 
linclude <nath.h> 
linclude "nui.h" 

void conv(float *, int, float int, float *) ; 
double sujDsq( float *, int) ; 

void Xb]nToZeee( float ^ibn, float *ieee, int n) 
/* 

This subroutine converts a 32 bit FP IBM number to a 32 bit I££E 
FP number 

The exponent of the IBM represents power of 16 against 2 for the 
ZSEE number and it is shifted one bit to the left. He also 
need to take care of the hidden bit of the IEEE fomat 

V 
{ 

register union { 
float f ; 
long 1; 

unsigned char c[43; 
) temp, tempi, temp2; 
unsigned char c; 
unsigned long swap test «■ 1; 
register int i; 

for (i=0; i<n; i++) { 

/* we do IBM to IEEE conversion */ 

tenp.f - ibaCi]; 

if (* (char *) 6swaptest} { 

/* byte swaping */ 

c « temp«c[0] ; 

temp.c[0] =■ temp.c[3]; 

temp,c(3] = c; 

c = temp . c [ 1 ] ; 

temp.c[l] = temp,c[23; 

temp, c [2) =• c; 

) 

if ((temp.l & 0x00800000) «=» 0) { /* test for high bit of IBM fraction 
temp2.1 « temp.l & 0x7f 000000? /* turn off sign bit 
tempi. 1 « temp.l & Oxff 000000; /* get exponent and sign 
ieee[i] « (tenp.f - tempi. f) * temp2.f * 0.125; 

) else ( 

temp2.1 = temp.l & -0x00800000; /* turn off high bit of IBM fraction 
tempi. 1 = temp.l & 0x7f 000000; /* get exponent and sign 
ieee[i] « temp2.f * tempi. f * 0.125; 
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IS 



45 



SO 



) 

) 



int flipit(int lnput_value) 
( 

union 
( 

int ival; 
char cval[2]; 
) a, b; 

a. ival » input_value; 

b. cvalJO] » a.cval[l]; 
b.cval[l] « a.cvaltO]; 
return (b. ival) ; 

} 



long longf lip (long input value) 
{ 

vmion 
( 

long Ival; 
char cval(4 3; 
20 ) a, b; 

a. Ival « input_value; 

b. cval[0] - a.cval[3]7 
b.cvalCl] « a.cvalC2]; 
b.cvalC23 » a.cvalfl]? 
b • cval [ 3 ] « a . cval [ 0 ) ; 

25 return (b. Ival) ; 



int normXcorr (float *x, int zpx, int nx, float *y, int zpy, int ny, float *c) 
{ 

20 int i, j? 

float ♦ptrx, *ptry, *ptrc, ♦lasty, *lastc? 
double factor; 

lasty » y + ny - 1; /* Set pointers. */ 

ptrx « x; 

35 for(i«0; i<nx; i++, ptrx++) /* Calculate cross-corr */ 

( /* values* V 

ptry a lasty; 

for(j=0; j<ny; ptry—) 
c[i+j] += (*ptrx) * (*ptry); 

40 
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) 

factor - sqrt(sumsq(x, nx) * suiasqCy, ny)); /* Calculate normali- 

/* zation factor. 

lastc = c + nx + ny - 2; 

for(ptrcac; ptrc<«lastc; ptrc++) /* Normalize cross-corr. 

♦ptrc /" factor; 

return(ny + zpx - zpy - 1); /* Return zero-lag index, 

) 

double sumsq (float *array, int num) 
( 

double squares » 0; 
int i; 

for{i«0; i<nuin; i++, array++) 

squares +« (*array) * (*array} ; 

return (squares) ; 

) 

void inax_and_lag (float *array, int numpts, float *max, int *lag) 
( 

float *ptr, *last; 
int pos; 

ptr = array; 

last » array + numpts; 

*nax - *ptr; 

pos e 0-( (nunpts-l)/2) ; 

*lag « pos; 

for(ptr++, POS++; ptr<last; ptr++, pos++) 
if(*ptr > *aax) 
{ 

*max " *ptr; 
*lag = pos; 

) 

> 

void shortXC (float *x, float *y, int npts, int nlags, float *cout) 
( 

double y_engy, *x_engy, *ptr_x_engy; 

float *factor, *ptr x, *ptr_y, *ptr c, *ptr fac, *last x, *last_fac; 
int i; 

x_engy - calloc(nlags, sizeof (double) ) ; /* Create arrays for trace 
factor « calloc(nlags, sizeof (float) ) ; /* segment energy and norm 

/* alization factors. 
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, /* Calculate energy of model 

y engy - Buiasq(y, npta) r trace segment, first trace 

:hV^r : !nl2trs^?^(U_.ngy * y.engy) ; /* segment, also first factpr 

/* calculate remaining trace 
ptr Cac - factor + 1; % energies and normaliztion 

ptr"x_engy - x_engy + l; factors. 
ptr_x - x; 

last X - X + npts; ^ 

fSrVS/fS-tLtJac^ Ptr.x_e„^*., ptr.x^, last...., 

* •ptr_x_engy- (MPtr_x_engy.l)) - j:?^!^^*/!^^!?^ , 
*ptr fac - (float) sqrt( (*ptr_x_engy) * y.engy) J 

free(x_engy) ; 

calculate correlation coefficients for nlags. 
?«t5-0. ptrWactor; l<nlags; i++, ptr_c++. ptr.fac**) 

^ ptr_x - X + i; 

ptr_y •= y; , 

last X « X + i + npts - i; 

for (7 ptr x<«last_x; ptr_x++. ptrj++) 

*ptr_c"+= (*Ptr.x) * (♦ptrjr); 
*ptr_c /= *ptr_fac; 

free ( factor) ; 

jirir.t^l^.^,V'^^^^^^^^ e<^al-len^n arrays. 

* float *ptrl, *ptr2, *last, sub=0, diff«0; 

ptrl « segl; 
ptr2 = seg2; 
last = ptrl + numpts; 

for(; ptrKlast; ptrl++, ptr2++) 

^ diff = *ptrl - *ptr2; 
if (diff<0) 

diff « 0 - diff; 
sun diff; 

} 

sum /= numpts; 
return (sum) ; 
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void conv(float *x, int nx, float *y, Int ny, float *c) 

^ float *ptrx, *ptry, *ptrc, •lastc; 
int i, j; 

/* Ensure that values in output array are zeros. */ 
ptrc «" c; 

lastc « c + nx + ny - 1; 
for(; ptrc<lastc; ptrc++) 
*ptrc - 0; 

ptrx - x; 

for(i-0> Knx; ptrx++) /* Calculate ♦/ 

( /* values. V 

ptry - y; 

for(j=0; j<ny; ptry++) 
cli+j] += (*ptrx) * (*ptry); 

} 

) 

void hilbt oper( float *oper, int npts) 
{ 

int i, nptsl; 

float rc, pi, *ptrl, *ptr2; 

pi « 4 * (float) atan(l.O); 
rc - 2 / pi; 
nptsl » (npts - 1) / 27 
ptrl = oper + nptsl +1; 
ptr2 = ptrl - 2; 

for(i=l; i<=nptsl; i+»2, ptrl+«2, ptr2-«2) 
( 

*ptrl « rc / i; 
♦ptr2 - -(*ptrl) ; 

) 



void phs_and_lag( float *real, int nreal, float *hbtoper, int noper, 
float *guad, float *phase, int *lag) 

^ float *ptr_real, *ptr_quad, *last, pea)c_env, envelope, convrt_deg; 
int pos; 

/* Convolve with Hilbert transfora operator to nake quadrature trace. */ 
conv(real, nreal, hbtoper, noper, quad); 

/* Set pointers, then calculate envelope and inst. phase at first sample. 
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ptr real - real; 

ptr~quad - quad + ( (noper - 1) / 2) ; 
last - real + nreal; 

convrt^deg « 180 / {4 * (float) atan(l.O)); 

pealc env = (float) sqrt ( ( *ptr.real * *ptr.real) + (^P^r ^uad * *ptr_quad) ) ; 
*phaie - ((float) atan2 (*ptr_quad, *ptr^real)) * convrt^^deg; 
pes • 0 - ((nreal - 1) / 2) ; 
*lag « pes; 

/* Search for peak of envelope, calculate phase at that point. */ 

for (ptr_real++ , ptr_,quad++ , pos++ ; ptr^reaKlast ; ptr_real-h+ , ptr_quad++ , pos++) 

^ envelope » (float) sqrt ( (♦ptr^real * *ptr_real) + (*ptr_quad * 
*ptr_quad) ) ; 
if (envelope > pea)c_env) 
{ 

peak_env « envelope; 
*lag « pos; 

★phase « ((float) atan2 (*ptr_guad, *ptr_real)) ♦ convrt_deg; 

) 

} 

) 

void correlate (float *x, int nx, float *y, int ny, float *c) 
{ 

int i, j; 

float *ptrx, *ptry, *ptrc, *lasty, *lastc; 

lasty = y + ny - l; /* Set pointers. */ 

ptrx - x; 

for(i-o; i<nx; i++, ptrx++) /* Calculate cross-corr V 

^ /* values. */ 

ptry = lasty; 

for(j-0; j<ny; ptry—) 
cCi+j] +- (*ptrx) * (*ptry); 

) 

> 



float dual_sign_bit( float *segl, float *seg2, int nunpts) 

^ float *ptrl, *ptr2, *last, seglslope, seg2slope; 
35 float disagree « 0; 

ptrl » segl +1; 
ptr2 « seg2 + 1; 
last " segl + nunpts; 

40 

for(; ptrKlast; ptrl++, ptr 2++) 

^ if(((*ptrl<0) && (*ptr2>0)) 11 {(*ptrl>0) fi* (*ptr2<0))) 
45 disagree++ ; 

seglslope - *ptrl - *(ptrl - 1); 

inaSLloS^f «*(Ieg2sloJi>0)) |l ( (seglslope>0) W (seg2slope<0) ) ) 
disagree++; 

) 

return (disagree) ; 

50 ) 
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SUBROUTINE F 

f include "nui.h" 

void q)cBort_asc( struct compare **, int, int) ; /* Ascending quicJcsort. */ 
void q)cBort_desc( struct compare **, int, int); /* Descending quicksort. */ 

void nuisort (struct compare **pointers, int numpts) 
( 

/* Sort elements: 

2 = xcorr coeff at zero lag 

3 m maximum xcorr coeff 

4 » average difference per sample 

5 " phase at max envelope of best x-corr 

6 =■ quality factor 

7 « # dual sign bit differences ♦/ 

75 gksort desc (pointers, ntmpts, 3); 

/* qksort asc (pointers, 15, 4); */ 

) 

void qksort_desc (struct compare **ptrs, int num, int offset) 
( /* Sort with a "quic)csort" algorithm. */ 

struct compare **1, **r, ***p, **i, *temp, **stack[50]; 

float v; 



20 



1 « ptrs; /* Initialize pointers to whole array, */ 

r " ptrs + num - 1; 
p « stack + 2; 

25 while (p >- stack +2) /* If there are partitions to do, and */ 

{ 

if(r > 1) /* if the partition is not zero-length, . */ 

( 

v » ♦((float *) (*r)+offset) ; /* Use last value in partition */ 

/* as "seed". */ 
i « 1; /* Set left pointer, */ 

j = r; /* Set right pointer. */ 



30 



do 
{ 

while(i < j && *((float *) (*i)+offset) >= v) 

i++; /* Increment pointers from left */ 

35 while(j > i && *((float *) (*j)+offset) <» v) 

j — ; /* and right. */ 

if(i < j) 
{ 

temp = *i; /* Swap out-of-order elements. */ 

*i = *j; 
*j « temp; 

) 

) while(i < j); 
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else 
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temp - *lj /* Put seed element in proper location */ 

*i . *r; /* in the array, separating two V 

*r »■ temp; /* partitions. V 

if((i - 1) > (r - i)) /* Determine which partition is larger. */ 
( 

*p - 1; /* V 

*(p+l) - 1 - 1; /* ^ V 

1 « i 1; /* Store starting and ending locations */ 

) /* of larger partition on stack 2md */ 

else /* reset left or right pointer to */ 

( /* proceed with shorter partition. */ 

*p » i + 1; /* */ 

*{p+l) - r; /* V 

r « i - 1; /♦ */ 

p 2; /* Increment stack pointer. V 



{ 

p -= 2; /* If done with short partition, */ 

1 a *p; /* retrieve earlier partition limits */ 

r ■ *(p+l) ; /* from stack for processing. */ 

) 



} 



void gksort_asc (struct compare **ptrs, int nun, int offset) 
( " /* Sort with a "quicksort" algorithm, 

struct compare **1, **r, ***p, **i, *temp, **stackC50]; 

25 float v; 

1 « ptrs; /* Initialize pointers to whole array. */ 

r «■ ptrs + num - 1; 
p - stack + 2; 



while (p >o stack + 2) /* If there are partitions to do, and */ 

if (r > 1) /* if the partition is not zero-length, */ 

{ 

V » *( (float *) (*r)+offset) ; /* Use last value in partition */ 

/* as "seed". */ 
i = 1; /♦ Set left pointer. */ 

35 j « r; /♦ Set right pointer. */ 

do 
{ 

while(i < j && *{(float *) (*i)+offset) <= v) 

i++; /* Increment pointers from left */ 

40 
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whilo(j > 1 tfc *((noat *) (*j)+offset) >» v) 
j.-; /* and right. */ 

if(i < j) 

^ temp « /* Swap out-of-order elements. */ 

♦1 - *j; 
*j * temp; 

) 

) while (i < j) ; 

temp = *i; /* P^t seed element in proper location */ 

*i B *r; /♦ in the array, separating two */ 

*r «= temp; /* partitions. */ 

if((i - 1) > (r - i)) /* Determine which partition is larger. */ 

^ *p - 1; /* V 

*(p+l) - i - 1; /* V 

l»i+l; /* store starting and ending locations */ 

) /* of larger partition on stack and */ 

else /* reset left or right pointer to */ 

{ /* proceed with shorter partition. */ 

*p - i + 1; /* V 

*(p+l) - r; /* V 

r = i - l; /* V 

p +s 2; /* Increment stack pointer. */ 

> 

else 

p -= 2; /* If done with short partition, */ 

1 « *p; /* retrieve earlier partition limits */ 

r « *(p+l); /* from stack for processing. */ 

} 
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* Get attribute data for display on graph. V 
fseelc(binfile, fileoftset, S££X_5ET) ; 

for{i«startindex, )c«0; i<endindex; i++, /* For each seis trace...*/ 

( 

if (seisattrib[i].error_flag — 1) 
{ 

fread(coDp, sizeof (struct conpare) , nimtnii, binf ile) ; 

continue; /♦ If top missing, skip all diplay. */ 

) 

/* Load "compare statistics for this trace and sort. */ 

for(j«0; j<nxantm; j++) /* build array of pointers.... */ 

comp_ptrs[j] ■ comp + j; 

fread(comp, sizeof (struct compare) , numtrm, binf ile) ; 

nuisort ( comp_ptrs , numtrm) ; 

/* Draw worst S points in green. */ 

hpos « grof f set + (PIXPERTRACE * k) ? 

setcolor(lO) ; 
for(J-5; j<10; j++) 

{ 

__setvieworg(PAYORG) ; 

attrib - iptattrib[comp_ptrs[ j ]->mod_cdp) .net_pay - minntpay; 
vpos - 0 - (int) (attrib * pixperpayunit) ; 
_rectangle(_GBORDER, hpos, vpos-1, hpos+1, vpos); 
_setvieworg(POKEORG) ; 

attrib - iptattrib[comp_ptrs( j ]->mod_cdp3 .net_por_ft - minntpore; 
vpos « 0 - (int) (attrib * pixperporeunit) ; 
_rectangle( GBORDER, hpoa, vpos-1, hpos+1, vpos); 
_setvieworgTHCPOREORG) ; 

attrib « iptattrib[comp_ptrs[ j ]->mod_cdp] ,hc_por_vol - minhcvol; 
vpos « 0 - (int) (attrib * pixperhcporeunit) ; 
_rectangle(_GBORDER, hpos, vpos-1, hpos+1, vpos); 

) 

Draw next four better points in yellow. */ 

_setcolor(14) ; 
for(j=l; j<S; j++) 
( 

_setvieworg(PAYORG) ; 

attrib = iptattrib[comp_ptrs( j ]->mod cdp].netjpay - minntpay; 
vpos - 0 - (int) (attrib * pixperpayunit); 
_rectangle(_GBORDER, hpos, vpos-1, hpos+1, vpos); 
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_Betvieworg(POREORG) ; 

attrib - iptattrib[conp_ptra[ j ]->inod_cdp] .net_por_ft - ninntpore 
vpos - 0 - (int) (attrib * pixperporeunit) ; *" 
_rectangle(_GBORDER, hpos, vpos-1, hpos+1, vpos); 
_Betvieworg(HCPORZORG) ; 

attrib « iptattr lb [ comp_ptr s [ j ] ->iaod_cdp } . hc_por_vol - ainhcvol ; 
vpos o 0 - (int) (attrib * pixperhcporeunit) ; 
_rectangle(_GBORDER, hpos, vpos-1, hpos+1, vpos); 



/♦ Draw best point in pin]c. */ 

8etcolor(12) ; 
^8etvievorg(PAY0RG) ; 

attrib - iptattr ib [ compjtrs [ 0 ] ->mod_cdp ] . netjay - ainntpay ; 
• vpos « 0 - (int) (attrib * pixperpayunlt) ; 
_rectangle(_GBORDER, hpos, vpos-1, hpos+1, vpos); 
26etvieworg(POREORG) ; 

attrib - iptattr ib [ comp_ptrs [ 0 ] ->mod_cdp ] . net __por_f t - ainntpore ; 
vpos » 0 - (int) (attrib * pixperporeunit) ; 
_rectangle( GBOROER, hpos, vpos*l, hpos+1, vpos); 
_setvieworgTHCPOREORG) ; 

attrib - iptattrib[conp_ptrs[0]->aod_cdp] .hc_por vol - ainhcvol; 
vpos « 0 - (int) (attrib * pixperhcporeunit) ; 
^rectangle (_GBORDER, hpos, vpos-l, hpos+1, vpos); 

* Plot averaged profiles through attributes. */ 

avgjpay - (3 * iptattrib[coap_ptrs[0]->nod_cdp] .netjay + 
2 * iptattribCcomp_ptrs(l]->ffiod_cdp] .net_pay + 
2 * iptattrib [ compjtrs [ 2 ] ->aod2cdp ] . netjay + 
2 * iptattrib [ comp^ptrs [ 3 ] ->nod_cdp ] . net_pay + 
iptattrib [ comp_ptrs { 4 ] ->aod_cdp ] . netjay + 
iptattrib [ coap^ptr s [ 5 ] ->aod_cdp ] . net^pay + 
iptattrib C coap_ptrs [ 6 ] ->aod_cdp ] . uetjay ) / 12 ; 

avg_pore - (3 * iptattribCcoap_ptrs[0]->aod_cdp] .net_j>or_ft + 
2 * iptattrib[coap_ptrs(l]->aod_cdp] .netjor^ft + 
2 * iptattribCcoap_ptrsC2 3->mod_cdp3 .netjor ft + 
2 * iptattrib[coap_ptrs(3]->aod_cdp] .netjor^ft + 
iptattribtcoap_ptrs[4]->aod_cdp] .netjor ft + 
iptattrib [ coap^tr s [ 5 ] ->BOd]]]cdp ] . netjor^ft + 
iptattrib [ coap_ptr s ( 6 ] ->aod_cdp ] . net_por_f t ) / 12 ; 

*vg_hcpore » (3 * iptattribCcoBp_ptrs[0]->aod_cdp3 •hc_por_yol + 
2 * iptattrib [coapj)trs[ 13 ->aod_cdp] .he jor_vol + 
2 * iptattribCcoap_ptrs[2]->aod2cdp] .hcjpor^vol + 
2 * iptattrib[coap_ptrs[3]->aod_cdp] .hc_j)or3vol + 
iptattrib [ coap^ptrs [ 4 ] ->aod_cdp 3 . he jpor_vol + 
iptattrib [coap_ptrs[ 53 ->aod2cdp 3 .hc_por_vol + 
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Iptattrlb [ comp_ptrs [ 6 ] ->nod_cdp ) • he jporjvol ) / 12 ; 



payvpoB - - (int) ( (avg_pay - ainntpay) * piacperpayunit) ; 

porovpos ■ - (int) ((avg_pore • minntpore) * pixperporeunit) ; 

hcporevpos - - (int) ((avg^hcpore - ninhcvol) * pixperhcporeunit) ; 

8etcolor(12) ; ^ 
If (i I- start index) 

{ 

_8etvieworg(PAyORG) ; 

_moveto(hpos - PIXPERTRACE, prevpayvpos) ; 
_lineto(hpos, payvpoa) ; 
_setvieworg(POREORG) ; 

jttovetp(hpos - PIXPERTRACE, prevporevpos) ; 

lineto(hpo8, porevpos) ; 
3Betvieworg(HCP0RE0RG) ; 

_moveto(hpos - PIXPERTRACE, prevhcvpos) ; . 
2lineto(hpo8, hcpor8vi>os) ; 

} 

prevpayvpos - payvpos; 
prevporevpos - porevpos; 
prevhcvpos ■ hcporevpos; 



Claims 

1. A method comprising the steps of: 

(1) providing at least one reference seismic trace corresponding to at least one reference lateral 
location and at least one nonreference seismic trace corresponding to at least one nonreference 
lateral location offset from said at least one reference lateral location, wherein each of the seismic 
traces results from the detection of the reflection of at least one seismic pulse as generated by at 
least one seismic source, and wherein each of the seismic traces includes a pair of reflection events 
respectively corresponding to the upper and lower boundaries of a subterranean layer of interest 

(2) providing a velocity log and a density log, which together comprise a reference log-pair, for the 
layer at said at least one reference lateral location so as to provide at least one reference-log pair 
having associated therewith at least one corresponding known value of at least one petrophysical 
property of the layer; 

(3) determining at least one reference reflection coefficient based on said at least one corresponding 
reference log-pair; 

(4) providing at least one seismic wavelet which is representative of said at least one seismic pulse 
at the layer and which when convolved with said at least one reference reflection coefficient 
produces at least one reference synthetic seismogram which approximates said at least one 
reference seismic trace; 

(5) deriving a set of a predetermined number of modified log-pairs, wherein each modified log-pair is 
different from one another and comprises a velocity log and a density log, each of which logs is a 
modified version of the respective velocity tog and density log of said at least one reference log-pair, 
and wherein each of the modified log-pairs has associated therewith at least one value of said at 
least one petrophysical property of the layer which is different than said at least one known value of 
said at least one petrophysical property and wherein such modified log-pairs are representative of 
probable variations of said at least one petrophysical property for said at least one nonreference 
lateral location; 

(6) determining a modified reflection coefficient based on the velocity log and density log of each of 
the modified log-pairs to thereby result in a number of modified reflection coefficients equivalent to 
the predetermined number of modified log-pairs; 

(7) convolving each of the modified reflection coefficients with said at least one seismic wavelet to 
produce a modified synthetic seismogram for each of the modified reflection coefficients to thereby 
result In a number of modified syntiietic seismograms equivalent to the predetermined number of 
modified log-pairs, each of the modified syntiietic seismograms having reflection events correspond- 
ing to the upper and lower boundaries of the layer and also having associated therewith at least one 
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value of said at least one petrophysical property associated with the corresponding modified log- 
pair; 

(8) comparing a connparison window of said at least one nonreference seismic trace to a conrespond- 
Ing comparison window of each of the modified synthetic seismograms, where each comparison 

5 window includes the reflection events corresponding to the upper and lower boundaries of the layer; 

(9) selecting those modified synthetic seismograms which match said at least one nonreference 
seismic trace sufficiently in step (8) to pass at least one predetermined matching threshold; 

(10) assigning to said at least one nonreference seismic trace the values of said at least one 
petrophysical property corresponding to the respective modified synthetic seismograms selected in 

70 step (9), thereby providing a range of values of said at least one petrophysical property associated 

with said at least one nonreference seismic trace and said at least one conresponding nonreference 
lateral location of the layer. 

2. A method as recited in claim 1 wherein said at least one petrophysical property is selected from the 
IS group consisting of layer thickness, porosity, lithological composition, water or hydrocarbon saturation, 

any derivative of one or more of the aforementioned properties, and combinations thereof. 

3. A method as recited in claim 2 wherein the velocity log and the density log of said at least one 
reference log-pair are obtained by deriving such logs from a reference set of logs comprising a 

20 lithological composition log, a porosity log, and a water or hydrocarbon saturation log which cor- 
responds to the reference lateral location, and wherein each of the modified log-pairs are derived from 
a modified set of logs of which at least one of such togs is a modified version of at least one of the logs 
of the reference set of logs. 

25 4, A method as recited in claim 3 wherein the velocity and density logs of said at least one reference log- 
pair and the velocity and density logs of the modified log-pairs are each a series of values expressible 
as a curve or function of time or depth. 

5. A method as recited in claim 4 wherein said at least one reference reflection coefficient and the 
30 modified reflection coefficients are each a series of values expressible as a curve or function of time or 

depth. 

6. A method as recited In claim 5 wherein said at least one seismic wavelet is shaped such that said at 
least one reference synthetic seismogram, resulting from convolution of said at least one seismic 

35 wavelet with said at least one reference reflection coefficient, matches said at least one reference 
seismic trace sufficiently to pass at least one predetermined matching threshold. 

7. A method as recited in claim 6 wherein in comparison step (8) the comparison employs crosscor- 
relation, 

40 

8. A method as recited in claim 6 wherein in comparison step (8) the comparison employs calculation of 
difference mismatch error. 

9. A method as recited in claim 6 wherein in comparison step (8) the comparison employs a combination 
45 of crosscorrelation and difference mismatch error. 

10. A method as recited in claim 6 wherein said at least one nonreference seismic trace comprises a 
plurality of nonreference seismic traces corresponding to a plurality of respective nonreference lateral 
locations. 

50 

11. A method as recited in one of the preceding claims comprising converting at least one of the 
seismograms obtained to a visible image thereof. 

1Z A method in accordance with one of the preceding claims comprising generating at least one seismic 
55 pulse by at least one seismic source and obtaining the reflection of said at least one seismic pulse and 
converting said reflection to said seismic trace. 
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